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' Temporal broadening of pulsar signals results from electron density fluctua- 

^ ■ tions in the interstellar medium that cause the radiation to travel along paths 

Tj- ■ of different lengths. The Gaussian theory of fluctuations predicts that the pulse 

^ ' temporal broadening should scale with the wavelength as A^, and with the dis- 

O ■ persion measure (corresponding to distance to the pulsar) as DM^. For large 

^ ■ dispersion measure, DM > 20pc/cm^, the observed scaling is \^DM^, contra- 

dieting the conventional theory. Although the problem has existed for 30 years, 
2 ' there has been no resolution to this paradox. 

^ ■ We suggest that scintillations for distant pulsars are caused by non-Gaussian, 

• spatially intermittent density fluctuations with a power-like probability distribu- 

tion. This probability distribution does not have a second moment in a large 
range of density fluctuations, and therefore the previously applied conventional 
Fokker-Planck theory does not hold. Instead, we propose to apply the theory of 
Levy distributions (so-called Levy flights). Using the scaling analysis (confirmed 
by numerical simulations of ray propagation) we show that the observed scaling 
is recovered for large DM, if the density differences, 5N, have Levy distribution 
decaying as \5N\~^/^. 

Subject headings: Turbulence — ISM: kinematics and dynamics 
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1. Introduction 



Intensity fluctuations of pulsars radiation are caused by scattering of radio waves by elec- 
tron density inhomogeneities in the interstellar medium. These fluctuations are a signature 
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of turbulent, non-equilibrium motion in the ISM, and as the phenomenon of turbulence itself, 
they have withstood full theoretical understanding for decades, see e.g. reviews by Sutton 
(1971) and Rickett (1977, 1990). Observationally, the presence of electron density fluctua- 
tions leads, among other effects, to temporal and angular broadening of the pulsar image. 
These two effects are naturally related - due to ffuctuations of the refraction index, different 
rays from a pulsar travel along paths of different shapes and the stronger the deviation of the 
path from the straight line, the broader the pulsar image and the larger the time broadening 
of the arriving signal. Denoting the angular width of the image as A^, and using sim- 
ple geometrical consideration, one estimates the arrival time broadening as Td ~ {AOyd/c, 
where d is the distance to the pulsar, and c is the speed of light, see more detailed discussion 
in (Blanford & Narayan 1985; Gwinn, Bartel, & Cordes 1993). 

A ray propagating through the interstellar medium encounters many randomly dis- 
tributed small "prisms" on its way, that make the scattering angle wander randomly. At 
each scattering event, the angle deflection is proportional to [see below], where A is the 
wavelength of the scattered radiation. Taking into account that the scattering angle is small 
and exhibits the standard Gaussian random walk, we estimate (A^)^ ~ A^rf, and the time 
delay scales as ~ X^d"^, where d plays the role of time in this random walk. The distance 
to the pulsar is proportional to the dispersion measure, DM, and therefore this relation 
can be checked experimentally. As has been consistently noted for more than 30 years, ob- 
served scahng of scintillations of distant pulsars, DM > 20pc/cm^, is far from this simple 
theoretical prediction, instead, it is well described by ~ X^DM^, see e.g., (Sutton 1971; 
Rickett 1977). Sutton proposed that scaling for longer hues of sight arose from dramatically 
increased probability of intersection with strongly scattering HII regions. In this sense, he 
proposed that rare, large events dominated the line-of-sight averages. 

The problem of scintillations was addressed by many authors who developed thor- 
ough analytical models, see the discussion in (Tatarskii & Zavorotnyi 1980; Rumsey 1975; 
Gochelashvily & Shishov 1975; Lee & Jokipii 1975a,b; Goodman & Narayan 1985; Blanford 
& Narayan 1985; Lithwick & Goldreich 2001). These models account for both smooth and 
non-smooth density fluctuations, the latter can arise from turbulent cascades. The main 
object of the theories is the so-called projected correlator of density fluctuations. Denote 
as N{y,t) the electron density, and 7V(x, t) = J (lzN{r,t) its projection perpendicular to 
the distance d. Here x is a two dimensional vector in the plane perpendicular to the line of 
sight, and 2; is a coordinate along the line of sight, i.e. r = (x, z). Note that these theories 
all assume that the distribution of projected density fluctuations is Gaussian. The density 
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and projected density correlators are related as: 

d d 

(iV(xi)iV(x2)) = ! [dzdz'{N{n)N{T2)), (1) 







where both fields inside the brackets are taken at the same time. Due to space homogeneity, 
these correlators depend only on the difference of the coordinates, e.g., (iV(ri)A^(r2)) = 
K{ri — r2). Assuming that the density fluctuations have flnite correlation length I, i.e. the k 
function decays fast for |ri — > we obtain 

oo 

(iV(xi)iV(x2)) = dj dzK{x^-X2,z) = dK(xi -X2). (2) 


It is easy to show that if in the inertial range of turbulent fluctuations, |r| -C /, the k function 
behaves as K{r) ki Nq[1 — B[r/l)°'], then the projected function is expanded as k{x,t) 
Nl[l - B{x/lY+% As an estimate, one has (A^^^ = ~ N^/l, and B and B are of 
the order 1. The analytical case corresponds to o; = 1, and in this case ~ X^dP. In a 
general case, the density held should not be analytic, and a^l. For example, Kolmogorov 
turbulence would imply a = 2/3. Such different possibilities have been exhaustively analyzed 
in the literature, see e.g. (Lcc & Jokipii 1975a,b; Goodman & Narayan 1985; Lambert & 
Rickctt 2000). Rigorous consideration shows that in the non-analytic case, the scaling of the 
broadening time changes. For a < 1, one obtains 

Td ~ X'^(.'^+3)/('^+^)d^<^+^y^'^+^\ (3) 

while for a more exotic case, a > 1, one gets 

Td~ A^/(^-")d(=^+")/(^-"). (4) 

In section 2 we present a simple derivation of these results. Since most observational data 
indicate that A-scaling is close to A"^, neither possibility provides enough freedom for changing 
the d-scaling from (P to d'^. 

In the present paper we propose a new model, that fully exploits the turbulent origin 
of the density fluctuations. We assume that the statistics of the density fluctuations is not 
Gaussian, but highly intermittent, and that the probability density function (PDF) of density 
differences has power-law decay, P{6N) ~ \5N\^^~^. If this power-law distribution does not 
have a second moment {(3 < 2), the Gaussian random walk approach does not work. Instead, 
we suggest to use the theory of Levy distributions, see (Shlesinger, Zaslavsky, & Frisch 
1995). Physically, the possibihty of power- law density distribution seems rather natural for 
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strong turbulent fluctuations. Indeed, the ISM turbulence can be near-sonic, i.e. velocity 
and density fields can develop shock discontinuities. Prom the theory of shock turbulence 
(Burgers turbulence) one knows that shocks or large negative velocity gradients have a 
power-law distribution, (Polyakov 1995; E et al 1997; Boldyrev 1998). Jump conditions on a 
shock then show that the velocity and density discontinuities arc proportional to each other, 
therefore density jumps may also have power-law distribution. Taking the Levy distribution 
of the density fluctuations as a working conjecture, we demonstrate that the scaling of the 
broadening time with respect to d is sensitive to the exponent of the distribution, /3, and the 
scaling ~ X^d^ is reproduced for /3 = 2/3. 

In the next section we review the ray-tracing model of pulse propagation, considered 
before by Williamson (1972, 1973); Blanford & Narayan (1985). In particular, we re-derive 
the results cited above for the Gaussian density fluctuations in a general, non-analj^ic case. 
In Section 3 we apply the model to the non-Gaussian, Levy distributed density fluctuations. 
We then numerically calculate the distribution of pulse-arrival times in the case of a smooth 
density field, and demonstrate that if P{SN) ~ \SN\~^/^, the width of this distribution 
changes with the distance to the pulsar as X'^d'^, in agreement with our scaling arguments. 
Conclusions and future research are outlined in Section 4. 



2. Ray-tracing method 

This method is applicable in the limit of geometrical optics, i.e. when the wave length 
is much smaller than the characteristic size of density inhomogeneities (Lifshitz, Landau, & 
Pitaevsky 1995). This rather effective method was applied to the problem of scintillations 
by Williamson (1972, 1973); Blanford & Narayan (1985); wc present it here in the form 
that allows us to apply it in the next section to Levy walks. In the limit considered, signal 
propagation can be characterized by rays, r(t), along which wave packets travel similar to 
particles obeying the following system of Hamilton equations: 

r = duj{k,r)/dk, 

k = —duj{k,r)/dr. (5) 

In this representation, u plays the role of Hamiltonian, uj^ = u!pe{r) + k'^c?, where oJp^ij') = 
47rN{r)e'^ /rrie is the electron plasma frequency and k is a wave vector. Differentiating the 
first equation in (5) with respect to t and using the second equation one obtains: 

= -27rc^X'^rodN{r)/dr, (6) 

where ro = e'^/irieC^ is the classical radius of electron. Taking into account that the ray 
propagates at small angles to the line of sight, chosen as the z-axis, we are interested in ray 
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displacement in the perpendicular, x direction, and instead of time we will use z variable, z — 
ct. Consider now two rays, separated by a vector Ax in the direction perpendicular to the 
2;-axis. As follows from (6), this vector obeys the following equation 

^ = Av. 

dz 

d(Av) ^ ^^9N^ 

dz as. 

where A — — 27rA^ro, and Av is an auhiliary variable having the meaning of velocity of 
beam spreading in the x direction, clearly A^ ~ Let us now assume that the electron 

density is a Gaussian random function with the correlation length /. Then, Av{z) is a 
Gaussian random walk, whose elementary time step has the length /. Since we are interested 
in very large propagation distances, z ^ I, and the scattering angles are very small, one can 
effectively assume that the random density is short-time correlated, i.e., the characteristic 
"z-time" of change of vectors Av and Ax is much larger than I. 

The diffusion coefficient for this random walk is 

and the diffusion is described by (A6')^ ~ Dz. We however observe that the diffusion 
coefficient depends of the distance Ax, and its behavior differs qualitatively for a < 1 and 
a > 1. In the first case, a < 1, diffusion is larger for smaller distances, therefore two rays 
are effectively attracting each other in the course of propagation. This means that at some 
point the geometrical ray picture will break down, and one needs to consider the effects of 
interference (interaction) of different rays. This happens when the beam is compressed to the 
size limited by the uncertainty condition in the perpendicular direction , kA6Ax ~ 1. Upon 
substituting A9 ~ D^^'^z^^'^, and using the expression for the diffusion coefficient (8), we can 
obtain the minimal size of contraction, and, equivalently, the diffraction angle corresponding 
to the aperture of this size. Assuming that the contraction happens at about half the distance 
between the pulsar and the Earth, z ~ d/2, we find: 

{Aef ^ [Ny,l-'^X'(^+'^d'] '^^"^'^ , a<l. (9) 

Recalling now that ~ [AO^d, we recover the result (3). In the second case, a > 1, the 
rays effectively repel, so geometrical optics does not break down. In this case Ax ~ Av z ~ 
£)i/22;3/2 'pjjjg equation gives 

{A9f - [NXl-^^X'^d^'^] '/^'""^ , 1 < a < 3, (10) 
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which agrees with the result (4). Both expressions give the same result for the analytic 
case, a = 1. The above standard results have been obtained by many authors and by a 
variety of different methods, see e.g., (Williamson 1972; Lee & Jokipii 1975a,b; Goodman 
& Narayan 1985; Blanford & Narayan 1985). As we mentioned in the introduction, neither 
one of the expressions (9) or (10) allows us to recover the observed scaling Td ~ A^rf^. In the 
next section we address the problem, assuming that the density- difference distribution has 
a slowly decaying power-law tail, such that the second moment of the distribution does not 
exist. In this case the diffusion approximation does not hold, and one needs to work directly 
with Eq. (7) to estabhsh the scahng of the probability of pulse arrival times. 



In previous section we implicitly used the central limit theorem, which states that the 
sum of many independent random variables has Gaussian distribution if second moments 
of these variables exist. More precisely, a convolution of many distribution functions that 
have second moments, converges to an appropriately rescaled Gaussian distribution. There- 
fore, the convolution of two Gaussian functions is a Gaussian function again. One can 
generalize this question for distribution functions without finite second moments: if their 
convolution converges, what is going to be the limit? The answer is the so-called Levy dis- 
tribution (Shlesinger, Zaslavsky, & Frisch 1995). As is the Gaussian distribution, the Levy 
distribution is stable: convolution of this distribution with itself gives the same distribution 
after proper rescaling. In other words, if two independent random variables are drawn from a 
Levy distribution, their sum has the same distribution, appropriately rescaled. Analogously 
to a Gaussian random walk, a sum of independent. Levy distributed random variables is 
called a Levy walk or Levy fiight. The latter name refiects the highly intermittent behavior 
of atypical Levy trajectory: it has sudden large jumps or "flights," see Fig. 1. Levy flights are 
common in completely different random systems and often replace diffusion in turbulent sys- 
tems. For example, a particle exhibiting a Brownian random motion in an equilibrium fluid, 
exhibits a Levy walk in a turbulent fluid. For a variety of further illustrations see (Shlesinger, 
Zaslavsky, & Frisch 1995). 

If a random variable y has a Levy probability density, P{y), then the Fourier transform 
of this distribution (the characteristic function) has the form: 



3. Levy model for scintillations 



oo 




(11) 



— oo 



where < (3 < 2, and C is some positive constant. For (3 — 2 we recover a Gaussian 
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distribution. This formula can be taken as the definition of a symmetric Levy walk. One 
can verify that P{y) ~ as \y\ oo. Of course, a distribution of a physical quantity 

usually has a second moment. This does not contradict our case, since the far tails of the 
PDF, which are not described by the Levy formula, make the dominant contribution to 
the second moment. However, if we are interested in effects caused by small fluctuations, 
y <^ yrms, it is the central part of the PDF that is important. 

The characteristic function of a convolution of n Levy distributions is just a product of 
n characteristic functions (11). We therefore conclude that the sum of n Levy distributed 
random variables has the distribution 



This is the demonstration of the convolution stability of the Levy distribution. Formula (12) 
teaches us that the displacement y of the Levy random walk scales with the number of steps 
as y ~ n^^^. In the Gaussian case, P — 2, we recover the well known diffusion result. 

We now would like to apply this result to our scintillation problem. Let us assume that 
the dimensionless density difference AN{x)/No has a Levy distribution with parameter (3. 
We then obtain from (5), Av ~ —AAN, and [compare this result to (8)!]: 



In this formula, a describes the scaling of the density fluctuations with distance, while (3 
is the exponent of the power-law decay of the density-difference probability distribution 
function. The scaling in Eq. (13) is understood not in the sense of averaging (the moments 
of A^^ of the order higher than f3 do not satisfy this scaling), but in the sense of scaling of the 
central part of the distribution Pz{A9). We now proceed exactly as we did in the derivation 
of formulae (9) and (10), and obtain for a < 1: 



(12) 




(13) 




l/(a+l) 



(14) 



and for 1 < q; < 3: 



r Ar4 4 ;2-2a-4//3 x 8 j2a-2+4//3] V (3"") 



(15) 



In the smooth (analytic) case, a — 1, the scaling of the time broadening is 



(16) 



We see that this scahng is sensitive to the exponent of the power distribution of the density 
fluctuations. This result was obtained by rather general arguments, and describes the scaling 
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of the arrival time distribution, rather than the moments of this distribution. Observations 
measure precisely the time width of the arriving signal, not its moments, i.e. they infer 
exactly the quantity corresponding to scaling (16). 

In the rest of this section we would like to verify the scaling (16), by numerical simulation 
of Eq. (7). Our simulations also provide the time-shape of the arriving signal. Let us assume 
that the distance to the pulsar, d, is much larger than the scale of an elementary scatter, 
i.e., n = d/l ^ 1, where n is the number of scattering events. At each scattering event, 
the angle of the ray changes hj 66 1, where 66 is a Levy-distributed random variable. 
[We denote 6N and 66 the characteristic changes of the density field and of the angle of 
propagation on one scattering segment of length / along the line of sight. This should not be 
confused with the changes of these variables between two different rays in a perpendicular 
plane x, denoted by A's.] The time delay (compared to the straight propagation) introduced 
by each scattering segment is Sra ~ W"^ /c. We need to find the probabihty distribution of 
the total travel time delay: 

2 

(17) 

assuming that each 66s [where 66s ~ ~A6N due to (7)] is distributed identically, indepen- 
dently, and according to the Levy law (11) with (3 = 2/3. This random variable can be 
generated in the following manner, see (Klaftcr, Zumofen, & Shlesinger 1995). Choose two 
positive numbers a > 6 > 1. Let 66 = 6qU with probability P{i) = (a — l)/(2a*"'"^), and 

66 = —6qU with the same probability, where i = 0, 1, 2 This is the so-called Weierstrass 

self-similar random walk, that can be considered as a discrete analog of the Levy walk with 
/3 = ln(a)/ln(&). 

In Fig. 2 we plot intensity of the arriving signal (the number of arriving rays) as a 
function of time. The arrival time was calculated with the aid of (17), where in the distri- 
bution of 66 we have chosen /3 = 2/3, 6 = 4, and = 0.0002. We considered the number 
of scattering events (the distance to the pulsar) to be ni — 100, n2 — 100 x 2^/^ 119, and 
^3 = 100 X 2~^/^ 84. Prom Fig. 2 one can see that the widths of the curves (estimated 
at the half of their maximum values) indeed differ by a factor of 2, as the scaling ~ A^d^ 
would predict for these distances. 

4. Conclusions 

In conclusion, we suggest a novel explanation for the observed scaling of time broadening 
of pulsar signals for large distances (large dispersion measures, DM > 20pc/cm^), ~ 



Td 



m=l 



m=l 



J2S6s 



s=l 
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A d . The central concept is that the density fluctuations in the intersteUar medium have a 
Levy probabihty distribution function that has power-law decay and does not have a second 
moment. The angle of pulse propagation, deviated by these density fluctuations, exhibits not 
a conventional Brownian motion, but rather a Levy flight. The exponent j3 is the parameter 
of the probability distribution of density differences, and the pulse broadening time is rather 
sensitive to it, as is described by our main formulae (14) and (15). The scaling ~ A^d^ 
is recovered for (3 = 2/3, i.e. for the \5N\~^/^ decay of the distribution function of density 
differences. This tail of the PDF appears as a result of turbulent density fragmentation, and 
it would be highly desirable to develop an analytical explanation for it. This is a concrete 
prediction of our model for the turbulence in the ISM, that can in principle be checked 
numerically. 
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Fig. 1. — A typical realization of a Levy random walk. The trajectory exhibits sudden 
large deviations, "flights." In the case of ray propagation through the ISM, the ray angle 
performs a Levy walk. Large angular deviations occur when the ray encounters regions 
of large electron density inhomogeneities, such as shocks or HII regions. (Angular scale is 
arbitrary.) 
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Fig. 2. — Numerical calculation of the number of arriving rays vs time (time units are 
arbitrary). We used formula (17), and the Levy distributed density fluctuations with f3 = 
2/3. We calculated arrival times of 10^ rays for three different distances to the source, 
Til = 100, n2 = 84 pa 100 x 2"^/^, and = 119 pa 100 x 2^/^. One observes that the width 
of the plot "n=84" is twice as small, and the width of the plot "n=119" is twice as large, as 
the width of the plot "n=100." This corresponds to the scaling ~ d^. 



